/*********************************************************************************
 * Software License Agreement (BSD License)
 * Copyright (c) 2012, Carlos D. Correa (info@ngraph.org)
 * All rights reserved.
 * 
 * Redistribution and use in source and binary forms, with or without
 * modification, are permitted provided that the following conditions are met:
 * 
 * 1. Redistributions of source code must retain the above copyright
 *    notice, this list of conditions and the following disclaimer.
 * 2. Redistributions in binary form must reproduce the above copyright
 *    notice, this list of conditions and the following disclaimer in the
 *    documentation and/or other materials provided with the distribution.
 * 3. Neither the name of netzen nor the
 *    names of its contributors may be used to endorse or promote products
 *    derived from this software without specific prior written permission.
 * 
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
 * ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
 * WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
 * DISCLAIMED. IN NO EVENT SHALL <COPYRIGHT HOLDER> BE LIABLE FOR ANY
 * DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
 * (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
 *  LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
 * ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
 * (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 *********************************************************************************/
/***************************************************************************
 *
 * Eigenvector Centrality Solver
 * based on IETL iterative solvers
 *
 * Copyright (C) 2009 by Carlos D. Correa <correac@cs.ucdavis.edu>
 *
 **************************************************************************/
#ifdef USINGLAPACK

#include "timing.h"

#include "node.h"
#include "centrality.h"


#include <boost/numeric/ublas/symmetric.hpp>
#include <boost/numeric/ublas/io.hpp>
#include <ietl/interface/ublas.h>
#include <ietl/power.h>
#include <ietl/vectorspace.h>
#include <ietl/rayleigh.h>
//#include <ietl/lanczos.h>
#include <ietl/iteration.h>
#include <boost/random.hpp>
#include <boost/limits.hpp>
//#include <ietl/solver.h>

typedef boost::numeric::ublas::matrix<double> Matrix; 
typedef boost::numeric::ublas::vector<double> Vector;
typedef ietl::vectorspace<Vector> Vecspace;
typedef boost::lagged_fibonacci607 Gen;

#include "lafnames.h"       
#include LA_GEN_MAT_DOUBLE_H
#include LA_VECTOR_DOUBLE_H 
#include LA_SYMM_MAT_DOUBLE_H
#include "blaspp.h"
#include LA_SOLVE_DOUBLE_H
#include LA_GENERATE_MAT_DOUBLE_H
#include LA_EXCEPTION_H
#include LA_UTIL_H
#include LA_SOLVE_DOUBLE_H 
#include <lasvd.h>

void My_Mat_Mat_Mult(int n, double *A, double *B, double *C, double a = 0) {
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++) {
      double sum = 0;
      for(int k=0;k<n;k++) {
	sum+=a*A[i*n+k]*B[k*n+j];
      }
      C[i*n+j] = sum;
    }
  }
}
 

inline void getDeriv(double *G, double *Qi, int i, int n) {
  for(int k=0;k<n*n;k++) Qi[k] = 0;

  for(int j=0;j<n;j++) {
    if(G[j*n+i]==0) continue;
    float sum = 0;
    for(int k=0;k<n;k++) {
      sum+=G[j*n+k];
    }
    if(sum>0) {
      for(int k=0;k<n;k++) {
	Qi[j*n+k] = -(G[j*n+k]/(sum+1.0) - G[j*n+k]/sum);
      }
      Qi[j*n+i] = -((G[j*n+i]+1.0)/(sum+1.0) - G[j*n+i]/sum);
    }
  }
}

void print(const char *title,  LaGenMatDouble &Y, int n) {
  printf("-------------     %s      --------------\n", title);
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++) {
      printf("%.4f\t", Y(i,j));
      //      printf("%.4f\t", Sig(i,j));
    }
    printf("\n");
  }
}

void print(const char *title,  LaVectorDouble &V, int n) {
  printf("-------------     %s      --------------\n", title);
  for(int i=0;i<n;i++) {
    printf("%.4f\n", V(i));
  }
}

MarkovCentrality::MarkovCentrality(Controller *controller):AnalysisTool("Centrality","MarkovCentrality",controller) {
    setPropertyBool("Derivatives", false);
    setPropertyStringList("Edge Weight", "Default", "edge");
    setPropertyString("Result", this->getName());	
}

void MarkovCentrality::compute() {
  assert(controller);
  // TO DO: set up these two via UI
  string resultProp = getPropertyString("Result");
  string resultDerivs = resultProp + ".Derivative";
	string resultDerivsT = resultProp + ".Derivative.Transpose";
	string resultDerivsM = resultProp + ".Derivative.Mutual";
	string resultDerivsB = resultProp + ".Derivative.Balance";
	string resultDerivsV = resultProp + ".Derivative.Magnitude";
	string resultDerivsN = resultProp + ".Derivative.Negative";
	
  Graph *g = controller->getGraph("graph");
  bool derivs = this->getPropertyBool("Derivatives");

  compute(g->nodes, derivs);
  controller->addDataProperty("node", resultProp, PROPERTY_FLOAT);
  g->computeNodeMinMaxFloat(resultProp);
  if(derivs) {
    controller->addDataProperty("edge",resultDerivs, PROPERTY_FLOAT);
    g->computeEdgeMinMaxFloat(resultDerivs);
  
	  controller->addDataProperty("edge",resultDerivsT, PROPERTY_FLOAT);
	  controller->addDataProperty("edge",resultDerivsM, PROPERTY_FLOAT);
	  controller->addDataProperty("edge",resultDerivsB, PROPERTY_FLOAT);
	  controller->addDataProperty("edge",resultDerivsV, PROPERTY_FLOAT);
	  controller->addDataProperty("edge",resultDerivsN, PROPERTY_FLOAT);

	  g->computeEdgeMinMaxFloat(resultDerivsT);
	  g->computeEdgeMinMaxFloat(resultDerivsM);
	  g->computeEdgeMinMaxFloat(resultDerivsB);
	  g->computeEdgeMinMaxFloat(resultDerivsV);
	  g->computeEdgeMinMaxFloat(resultDerivsN);
  }

  /*
    controller->exportFloat(graph,"markov.float", resultProp);
    if(derivs) {
    controller->exportImplicitEdgeFloat(graph, "markovDerivs.float", resultDerivs);
    graph.computeImplicitEdgeMinMaxFloat(resultDerivs);
    }
  */
}


int MarkovCentrality::compute(vector<Node*> &nodes, bool needDerivs) {
  int n = nodes.size();
  double *G = new double[n*n];
  double *data = new double[n*n];
  double *characteristic = new double[n*n];
  double *qi = new double[n*n];

  string property = getPropertyString("Edge Weight");
  string resultProp = getPropertyString("Result");
  string resultDerivs = resultProp + ".Derivative";
	string resultDerivsT = resultProp + ".Derivative.Transpose";
	string resultDerivsM = resultProp + ".Derivative.Mutual";
	string resultDerivsB = resultProp + ".Derivative.Balance";
	string resultDerivsV = resultProp + ".Derivative.Magnitude";
	string resultDerivsN = resultProp + ".Derivative.Negative";
	
  if(listener) listener->setProgress(0);

  for(int k=0;k<n*n;k++) {G[k] = 0; data[k] = 0;}

  for(int j=0;j<n;j++) {
    for(unsigned int k=0;k<nodes[j]->edges.size();k++) {
      assert(nodes[j]->edges[k]);
      assert(nodes[j]->edges[k]->getNode1());
      int i = nodes[j]->edges[k]->getNode1()->getIndex();
      float weight = 1.0;
      if(nodes[j]->edges[k]->hasProperty(property)) {
	weight = nodes[j]->edges[k]->getPropertyFloat(property);
      }
      G[j*n + i] += weight;
    }
  }

  for(int j=0;j<n;j++) {
    float sum = 0;
    for(int i=0;i<n;i++) {
      sum += G[j*n+i];
    }
    if(sum>0) {
      for(int i=0;i<n;i++) {
	data[j*n+i] = G[j*n+i]/sum;
      }
    }
  }

  for(int j=0;j<n;j++) {
    for(int i=0;i<n;i++) {
      double identity = (i==j)? 1.0:0.0;
      characteristic[j*n+i] = identity - data[j*n+i];
    }
  }


  timestamp t1 = now();

  // Initialize Matrix
  Matrix A(n,n);
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++) {
      A(i,j) = data[j*n+i];
    }
  }

  ietl::vectorspace<Vector> vec(n);
  boost::lagged_fibonacci607 gen;  
  int max_iter = 200;
  double rel_tol = 50.*std::numeric_limits<double>::epsilon();
  double abs_tol = std::numeric_limits<double>::epsilon();

  //double rel_tol = 1.*std::numeric_limits<double>::epsilon();
  //double abs_tol = std::numeric_limits<double>::epsilon();
  ietl::basic_iteration<double> iter(max_iter, rel_tol, abs_tol);
  
  print("Computing markov Power method max_iters = %d rel_tol = %g\n", max_iter, rel_tol);
  std::pair<double,Vector> result = ietl::power(A, gen, iter, vec); 

  Vector error = ietl::new_vector(vec);
  ietl::mult(A,result.second,error);
  error -= result.first*result.second;
  std::cout << "error: " << ietl::two_norm(error) << std::endl;

  print("Eigenvalue = %f\n", result.first);
  LaVectorDouble X(n);
  for(int i=0;i<n;i++) {
    X(i) = result.second[i]; 
    //printf("%g\n", X(i));
  }
  
  LaGenMatDouble Y(n,n);
  LaGenMatDouble Z(n,n);

  LaGenMatDouble VT(n,n);
  LaGenMatDouble U(n,n);
  LaGenMatDouble T(n,n);
  LaVectorDouble S(n);
  LaGenMatDouble Sig(n,n);
     
  double epsilon = 1e-10;
  print("Z = Inverse of Y\n");
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++) {
      float I = (i==j)? 1.0:0.0;
      Y(i,j) = I - A(j,i) - X(j);
    }
  }


  LaSVD_IP(Y, S, U, VT);
  for(int i=0;i<n;i++) {
    for(int j=0;j<n;j++) Sig(i,j) = 0.0;
    Sig(i,i) = fabs(S(i))>epsilon? 1/S(i): 0;
  }
  Blas_Mat_Trans_Mat_Mult(VT, Sig, T);
  Blas_Mat_Mat_Trans_Mult(T, U, Z);
    

  LaGenMatDouble IZE(n,n);
  LaGenMatDouble D(n,n);
  for(int i=0;i<n;i++) {
    D(i,i) = (1.0/X(i) + 0);
    for(int j=0;j<n;j++) {
      if(i!=j) D(i,j) = 0;
      double I = (i==j)? 1.0:0.0;
      IZE(i,j) = I - Z(i,j) + Z(j,j);
    }
  }


  LaGenMatDouble M(n,n);
  Blas_Mat_Mat_Mult(IZE, D, M);


  //print("D", D, n);


  LaVectorDouble Imp(n);
  for(int j=0;j<n;j++) {
    double sum = 0;
    for(int i=0;i<n;i++) {
      sum+=M(i,j);
    }
    Imp(j) = (float)n/sum;
  }

  timestamp t2 = now();
  

  FILE *fp1 = fopen("imp.txt", "w");
  for(int i=0;i<n;i++) {
    nodes[i]->setPropertyFloat(resultProp, Imp(i));
    nodes[i]->createImplicitEdges(n, nodes);
    fprintf(fp1, "%g\n", Imp(i));
  }
  fclose(fp1);


  if(needDerivs) {
    LaGenMatDouble Q(characteristic, n, n);
    LaGenMatDouble Qinv(n, n);
    LaGenMatDouble U(n,n);
    LaGenMatDouble VT(n,n);
    LaGenMatDouble Sig(n,n);
    LaGenMatDouble T(n,n);
    LaVectorDouble S(n);
    print("PseudoInverse of Q\n");
    LaSVD_IP(Q, S, U, VT);
    for(int i=0;i<n;i++) {
		for(int j=0;j<n;j++) {
			Sig(i,j) = 0;
		}
		Sig(i,i) = fabs(S(i))>1e-6? 1/S(i): 0;
    }
    Blas_Mat_Trans_Mat_Mult(VT, Sig, T);
    Blas_Mat_Mat_Trans_Mult(T, U, Qinv);
   
    //    print("Qinv", Qinv, n);

    LaVectorDouble QiX(n);
    LaVectorDouble dPi(n);
    LaGenMatDouble Qed(n,n);
    LaGenMatDouble Zi(n,n);
    LaGenMatDouble T1(n,n);
    LaGenMatDouble Di(n,n);
    LaGenMatDouble Mi1(n,n);
    LaGenMatDouble Mi2(n,n);
    LaGenMatDouble Mi(n,n);
    LaGenMatDouble ZiEZdg(n,n);
    LaVectorDouble Derivi(n);

    for(int k=0;k<n;k++) {
      if(listener) listener->setProgress((float) k/(float) n);
      print("Computing derivatives %d(%d)\r", k, n); fflush(stdout);
      getDeriv(G, qi, k, n);
      LaGenMatDouble Qi(qi, n, n);

	  Blas_Mat_Vec_Mult(Qi, X, QiX, 1.0);
      Blas_Mat_Vec_Mult(Qinv, QiX, dPi, -1.0);

      for(int i=0;i<n;i++) {
		for(int j=0;j<n;j++) {
		  Qed(i,j) = Qi(j,i) - dPi(j);
		}
      }
      
      Blas_Mat_Mat_Mult(Z, Qed, T1, -1.0);
      Blas_Mat_Mat_Mult(T1, Z, Zi);

	  for(int i=0;i<n;i++) {
		Di(i,i) = -1.0/(X(i)*X(i)) * dPi(i);
		for(int j=0;j<n;j++) {
		  if(i!=j) Di(i,j) = 0;
		}
      }
      for(int i=0;i<n;i++) {
		for(int j=0;j<n;j++) {
		  Mi(i,j) = Di(j,j)* IZE(i,j) +   D(j,j) * ( -Zi(i,j) + Zi(j,j) );
		}
      }
      for(int j=0;j<n;j++) {
		double sum = 0;
		for(int i=0;i<n;i++) {
		  sum+=Mi(i,j);
		}
		double avg = sum/(double) n;
		Derivi(j) = -(Imp(j)*Imp(j))*avg;
      }
      for(int j=0;j<n;j++) {
		  nodes[k]->getImplicitEdge(j)->setPropertySigned(resultDerivs, Derivi(j));
		  if(nodes[k]->hasEdge(j)) nodes[k]->getEdge(j)->setPropertySigned(resultDerivs, Derivi(j));		

		  nodes[j]->getImplicitEdge(k)->setPropertySigned(resultDerivsT, Derivi(j));
		  if(nodes[j]->hasEdge(k)) nodes[j]->getEdge(k)->setPropertySigned(resultDerivsT, Derivi(j));
		  
		  nodes[k]->getImplicitEdge(j)->setPropertySigned(resultDerivsV, fabs(Derivi(j)));
		  if(nodes[k]->hasEdge(j)) nodes[k]->getEdge(j)->setPropertySigned(resultDerivsV, fabs(Derivi(j)));
		  
		  nodes[k]->getImplicitEdge(j)->setPropertySigned(resultDerivsN, -Derivi(j));
		  if(nodes[k]->hasEdge(j)) nodes[k]->getEdge(j)->setPropertySigned(resultDerivsN, -Derivi(j));
      }
	}
	  
	for(int i=0;i<n;i++) {
	  for(int j=i;j<n;j++) {
		  float di, dj;
		  di = nodes[i]->getImplicitEdge(j)->getPropertyFloat(resultDerivs);
		  dj = nodes[j]->getImplicitEdge(i)->getPropertyFloat(resultDerivs);
		  float mutual = fabs(di)>fabs(dj)? di: dj; 
		  float balance = fabs(di-dj);
		  
		  nodes[i]->getImplicitEdge(j)->setPropertySigned(resultDerivsM, mutual);
		  nodes[j]->getImplicitEdge(i)->setPropertySigned(resultDerivsM, mutual);
		  
		  nodes[i]->getImplicitEdge(j)->setPropertySigned(resultDerivsB, balance);
		  nodes[j]->getImplicitEdge(i)->setPropertySigned(resultDerivsB, balance);
		  
		  if(nodes[i]->hasEdge(j)) nodes[i]->getEdge(j)->setPropertySigned(resultDerivsM, mutual);
		  if(nodes[j]->hasEdge(i)) nodes[j]->getEdge(i)->setPropertySigned(resultDerivsM, mutual);
		  
		  if(nodes[i]->hasEdge(j)) nodes[i]->getEdge(j)->setPropertySigned(resultDerivsB, balance);
		  if(nodes[j]->hasEdge(i)) nodes[j]->getEdge(i)->setPropertySigned(resultDerivsB, balance);
		  
	  }
	}	  
  
  }
  printf("Computing derivatives DONE\n");
      
  timestamp t3 = now();
  

  
 
  delete G;
  delete data;
  delete characteristic;
  delete qi;
  printf("[MARKOV CENTRALITY] Time  %f sec\n", t2 - t1);
  printf("[MARKOV DERIVS] Time  %f sec\n", t3 - t2);
  if(listener) listener->done();
  return 1;
}

#endif// USINGLAPACK


